Hexagons and Interfaces in a Vibrated Granular Layer 
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The order parameter model based on parametric Ginzburg-Landau equation is used to describe 
high acceleration patterns in vibrated layer of granular material. At large amplitude of driving 
both hexagons and interfaces emerge. Transverse instability leading to formation of "decorated" 
interfaces and labyrinthine patterns, is found. Additional sub-harmonic forcing leads to controlled 
interface motion. 
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Driven granular materials often manifest collective 
fluid-like behavior: convection, surface waves, and pat- 
tern formation (see e.g. jD). One of the most fascinat- 
ing examples of this collective dynamics is the appear- 
ance of long-range coherent patterns and localized exci- 
tations in vertically- vibrated thin granular layers 
The particular pattern is determined by the interplay 
between driving frequency / and acceleration of the con- 
tainer r = 4ir 2 Af 2 / g [A is the amplitude of oscillations, 
g is the gravity acceleration) (2j||]. Patterns appear at 
r rj 2.4 almost independently of the driving frequency /. 
At small frequencies / < /* || the transition is subcrit- 
ical (hysteretic), leading to formation of either squares 
or localized excitations such as individual oscillons and 
various bound states of oscillons. For higher frequencies 
/ > /* the transition is supercritical, and stripes appear 
first. Both squares and stripes oscillate at half of the 
driving frequency f/2. At higher acceleration (r > 4), 
stripes and squares become unstable, and hexagons ap- 
pear instead. Further increase of acceleration at T » 4.5 
converts hexagons into a domain-like structure of flat lay- 
ers oscillating with frequency f/2 with opposite phases. 
Depending on parameters, interfaces separating flat do- 
mains, are either smooth or "decorated" by periodic un- 
dulations. Undulations slowly grow with time and even- 
tually become quasi-periodic labyrinthine patterns. For 
r > 5.7 various quarter-harmonic patters emerge. 

The problem of pattern formation in thin layers of 
granular material was studied by several groups. While 
molecular dynamics simulations |J and hydrodynamic 
and phenomenological models [fjj reproduced certain ex- 
perimental features, neither of them offered a systematic 
description of the whole rich variety of the observed phe- 
nomena. In Ref. ||] we introduced the order parameter 
characterizing the complex amplitude of sub-harmonic 
oscillations. The equations of motion following from the 
symmetry arguments and mass conservation reproduced 
essential phenomenology of patterns near the threshold 
of primary bifurcation: stripes, squares and oscillons. 

In this Letter we describe high acceleration patterns 
on the basis the order parameter model. We show that 



at large amplitude of driving both hexagons and inter- 
faces emerge. We find morphological instability leading 
to formation of "decorated" interfaces and labyrinthine 
patterns. We develop the description of the labyrinths 
in terms of nonlocal contour dynamics. Additional sub- 
harmonic forcing leads to the motion of the interface with 
the direction controlled by relative phases of harmonic 
and sub-harmonic components of forcing. 

The essential of the model is the order parameter 
equation for the complex amplitude ip of parametric layer 
oscillations h = ip exp(j7r/i) + c.c. at the frequency f/2. 
In general, ip is coupled with the the average thickness 
of the layer p. However, at high frequencies (/ > /*) 
the coupling between p and ip becomes unimportant, and 
the model can be reduced to a single order-parameter 
equation 



d t i> = yip* - (1 - iw)ip + (1 + ib)V 2 ip 



2 V> (1) 



Eq. (|lj) describes the evolution of the order parameter for 
the parametric instability in spatially-extended systems 
(see 0,0). Linear terms in this equation are obtained 
from the complex growth rate for infinitesimal periodic 
layer perturbations h ~ exp[A(fc)i + ikx]. Expanding 
A(fc) for small k, and keeping only two leading terms 
in the expansion A(k) = —Aq — Aifc 2 we reproduce the 
linear terms in Eq. (|l|), where b — lmA\j ' ReA\ and 
lj = (f2o — 7r /)/^eAo, where fio = —ImA^. The term 
jip* characterizes the effect of driving at the resonance 
frequency. The term — | , 0| 2 V ; accounts for nonlinear sat- 
uration of waves at finite amplitude. 

It is convenient to shift the phase of the complex order 
parameter via ip = ipexp(i<f>) with sin2</> = w/'y. The 
equations for real and imaginary part ip = A + iB are: 

d t A = (s - l)A - 2loB - (A 2 + B 2 )A + V 2 (A - bB) (2) 



d t B = -(s + \)B - {A 2 + B 2 )B + V 2 {B + bA) 



(3) 



where s 2 = 7 2 — to 2 . At s < 1, Eqs. (0),(@) has only one 



trivial uniform state A = 0, B 
uniform states appear, A = Aq 



0, At s > 1, two new 
0. 



±Vs^T, B = 0. The 
onset of these states corresponds to the period doubling 
of the layer flights sequence, observed in experiments ^ 
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and predicted by the simple inelastic ball model [p|,|l2"|]. 
Signs ± reflect two relative phases of layer flights with 
respect to container vibrations. 

First we analyze the stability of the state A = Aq , B = 
with respect to perturbations with wavenumber k, 
(A, B) = (A ,0) + (U k ,V k ) exp[\(k)t + ikx\. The uniform 
state loses its stability with respect to periodic modula- 
tions with the wavenumber k c at s < s c (correspondingly, 
7 < 7 C ), where 



y/(l + LU 2 )(l + b 2 ) - U + b 



2b 

2s - 1 - ub 
1 + b 2 



(4) 
(5) 



Small perturbations with all directions of the wavevector 
grow with the same rate. The resultant selected pat- 
tern is determined by the nonlinear competition between 
the modes. In the presence of the reflection symmetry 
tp — > —ip, quadratic nonlinearity is absent, and cubic non- 
linearity favors stripes corresponding to a single mode. 
Near the fixed points A = Ao , B = the reflection sym- 
metry for perturbations U — > — U, V ^ — V is broken, 
and hexagons emerge at the threshold of instability. To 
clarify this point we perform weakly-nonlinear analysis 
of Eqs. (@),(||) for s = s c - e, and e < 1. At e — > 0, the 
variables U and V are related as in linear system: 

(U, V) = (1, V )*,V= [2(s c - 1) + k 2 c ] /(bk 2 c - 2w) 
W = Aj exp[zkr] + c.c, |fc[ = fc e (6) 

Corresponding adjoint eigenvector is 

(U+,V+) = (1,77+) , 77+ = - [2(s B 1) + k 2 } /bk 2 (7) 

Substituting Eq. (j6|) into Eqs. (|^),(^) and performing 
the orthogonalization, we obtain equations for the slowly- 
varying complex amplitudes Aj, j = 1,2,3 (we assume 
only three waves with triangular symmetry, favored by 
quadratic nonlinearity) 



d t Aj = 2eAj + a 2 A) +1 A*_ x 

-a 3 (\AA 2 + 2(\A^ 1 \ 2 + \A J+1 \ 2 ))A J 

where the coefficients 02,03 are 



a 2 = 2A 2 



1 + rjrj+ 



, a 3 = 3(1 + T]T] + ) 



(8) 



(9) 



Eqs. (g) are well-studied (see M]). There are three 
critical values of e: ca — — a\/AQa 3 , — a 2 /2a 3l and 
eg = 2a 2 /a 3 . The hexagons are stable for < e < eg, 
and the stripes are stable for e > e^. Thus, near s = s c 
the model exhibit stable hexagons Ej. Since we have 
two symmetric fixed points, both up- and down-hexagons 
co-exist. For smaller s stripes are stable, and for larger 
s, flat layers are stable, in agreement with observations 



H|16|]. The above analysis requires the values of ca,b,r to 
be small. For parameters u>,b = O(l), this requirement 
is satisfied for ea, but not for 63,11- The estimates can 
be improved by substituting Ao at s = s c + e instead of 
Aq(s c ) in (^|). The resulting range of stable hexagons is 
plotted in the phase diagram (a;, 7) (see Fig. |l]). 

At s > 1, Eqs. (||),(H) have an interface solution 
connecting two uniform states Aq — ±ys— 1, B = 
0. For b = this solution is of the form A = 
Aq tanh(Aoa;/2), B = 0. For b ^ the solution is avail- 
able only numerically. In order to investigate the stabil- 
ity of the interface, we consider the perturbed solution 
(A, B) = (A ,B ) + (A(x),B(x))exp[X(k)t + iky}. For 
A, B we obtain linear equations 



bk' 



B 

-A 



(10) 



s - 1 - ZAl - Bl + <9 2 , -2w - 2A B - bd 2 

-2A B + bd 2 , -s - 1 - A% - 3B 2 + 9 2 



In order to determine the spectrum of eigenvalues X(k) 
one has to solve Eq. (|l^) along with stationary Eqs. 
(||),||). Using numerical matching-shooting technique, 
we have obtained positive eigenvalues corresponding to 
interface instability. This instability is confirmed by di- 
rect numerical simulations of Eq. (0) . An example of the 
evolution of slightly perturbed interface is shown in Fig. 
||. Small perturbations grow to form a "decorated" in- 
terface similar to [0 , with time these decorations evolve 
slowly and eventually form a labyrinthine pattern [fl6f . 

The neutral curve for this instability can be deter- 
mined as follows. Numerical analysis shows that at the 
threshold the most unstable wavenumber is k — and 
we can expect that for k — > OA ~ k 2 . Expanding 
Eqs.© in power series of k 2 : (A, B) = (A^,B^) + 
k 2 {A^>,B^) + .... in the zeroth order in k we obtain 
L(A^°\B^°) = 0. The corresponding solution is the 
translation mode ,4 (0) = d x A (x), B^ = d x B (x). In 
the first order in k 2 we arrive at the linear inhomogc- 
neous problem 



(Hk) + k 2 )(^Z) +bk2 



B (o) 
-A(°) 



(11) 



A bounded solution to Eq. ( |11| ) exists if the r.h.s. is 
orthogonal to the localized mode of the adjoint operator 
A + ,B + . The orthogonality condition fixes the relation 
between A and k: 



X = -(l + (3)k 2 ,(3 



bJZo( AmB+ ~B^A+)dx 
(A(°)A+ + B(°)B+)dx 



(12) 



The instability, corresponding to the negative surface 
tension of the interface, onsets at = — 1. The neu- 
tral curve is shown in Fig. [I]. This instability leads to so 



2 



called "decorated" interfaces (see Fig. 0,a). At the non- 
linear stage the undulations grow and form labyrinthine 
patterns (Fig. g, b-d). The evolution of the inter- 
face can be investigated near the line s = 1 (see Fig. 
|l|). In the vicinity of this line A ~ (s — l) 1 / 2 and 
B ~ (s — l) 3 / 2 <C A. in the leading order we can ob- 
tain from Eq. (g) B = 6V 2 A/2, and Eq. (g) yields 

fr 2 

<M = (s - 1)A - A 3 + (1 - w&)V 2 ,4 - — V 4 A (13) 

Rescaling the variables t — > (s — l)t : A — > (s — l) -1 / 2 ^, 
.x -> (2(s— I)/?? 2 ) 1 / 4 , we can reduce Eq. (|l|) to the Swift- 
Hohcnberg equation (SHE) 



d t A = A - A d - SV 2 A - V 4 A 



(14) 



where S = (uib— l)y / 2/((s — 1)6 2 ). This equation is sim- 
pler than Eq. ([|) , however it captures many essential fea- 
tures of the original system dynamics, including existence 
and stability of stripes and hexagons in different parame- 
ter regions (see flfi}] ), existence of the interface solutions, 
interface instability and emergence of labyrinthine pat- 
terns. Indeed, simple analysis shows that the growth 
rate of the instability of the uniform state A = 1 as a 
function of the perturbation wavenumber is determined 
by the formula A(fc) = — 2 + 8k 2 — fc 4 , so it becomes un- 
stable at S > S c — 2\[2 at critical wavenumber k c = \/2. 
As in the original model, near the threshold of this insta- 
bility, subcritical hexagonal patterns are preferred. In- 
terface stability can also be analyzed more simply as the 
linearized operator corresponding to the model ( |i"4| ) is 
self-adjoint. The threshold value 6^ = 1.011 is obtained 
from the following solvability condition 



( 5 th A 0x 



2Ai 



dx = 



(15) 



We derived the equation for the position of the curved 
interface R parameterized by its arclength £. The solu- 
tion is sought for in the form A = Aq(ii ■ (r — R(£))) + 
W(r), where n is the unit normal vector, Aq is the ID 
interface solution and W is the correction due to interac- 
tion with the remote parts of the interface. This correc- 
tion can be obtained from the linearized SHE in the far 
field. Substituting this solution in the SHE (Q), we ob- 
tained the equation for the interface velocity c n — n<9 t R, 

c n = -™ + ^^dOTi(x|R(0-R(£)l) + cc. (16) 

where k is the local interface curvature, p — J A 2 dx 
is an interface "mass" per unit length, and a = 
—5 + 2 J A 2 x dx/fi is the surface tension, £ = (R(£) — 

RO/|R(0-R(OxR C (0 and X 2 = 6/2-y/P/4=2. 
Coefficient C is determined from the matching W(x) 
with the asymptotic of one-dimensional interface solu- 
tion. Eq. (Il6|) is similar to that of Ref. [18[ . The first term 



in the r.h.s. of ( [Tq ) at a > describes the destabilizing 
effect of negative surface tension, and the Biot-Savart 
type integral accounts for non-local interactions among 
parts of the interface. According to Ref. combina- 
tion of these effects gives rise to labyrinthine patterns. 
Since \ here is complex, function Kg oscillates, and in- 
terfaces form bound states at the distance I ~ 0(x _1 ). 
The characteristic wavelength of labyrintine patterns is 
I and is different from the most unstable wavelength of 
the original interface instability. This difference is evi- 
dent from Fig. g which shows snapshots of the interface 
dynamics. 

In the region of stability (/3 > —1) for the original 
model Eq. (|l|) the interface is stationary due to symme- 
try. However, if the plate oscillates with two frequencies, 
/ and //2, the symmetry between two states connected 
by the interface, is broken, and interface moves. The ve- 
locity of interface motion depends on the relative phase 
of the sub-harmonic forcing with respect to the forcing 
at /. This effect can be described by the additional term 
qe l ® in Eq. (|j), where q characterizes the amplitude of 
the sub-harmonic pumping, and $ determines its rela- 
tive phase. For small q, we look for moving interface 
solution in the form ip = ipo(x — vt) + qipi(x — vt) + .... 
and v = 0(q). Solvability condition yields the following 
expression for the interface velocity 

cos$ f A+dx + sin$ f B + dx , . 

v = — a (17) 

q J(A+d x A + B+d x B )dx K ' 

The explicit answer is possible to obtain for b — when 
A + = d x Ao, and <& = 0, tt, which yields the interface 



velocity v = +~^qA Q 2 — T^q(s — 1) 4 . In general, A, B, 



A + , B + , and hence v, can be found numerically. The 
interface velocity as function of q, $ is shown in Fig. ||. 

We have shown that large acceleration patterns are 
captured by the generic parametric Ginzburg-Landau 
equation (|l|) . The structure of the phase diagram Fig.|l|is 
qualitatively similar to that of experiments Refs. |p|, |l6|jl7| 
for high frequencies of vibration. Increasing vibra- 
tion amplitude leads to transition from a trivial state 
to stripes, hexagons, decorated interfaces, and finally, 
to stable interfaces. In experiments, at yet higher T, 
quarter-harmonic patterns appear, however these pat- 
terns are not described by our model. Transition from 
unstable to stable interfaces also occurs with decreas- 
ing oj (increasing vibration frequency /), in agreement 
with Ref. 0,0. In our model, the interface instability 
leads to labyrinthine patterns, however it seems that in 
experiments this instability sometimes saturates to pro- 
vide stationary "decorations" . Presumably, this satura- 
tion is caused by the weak interaction with the density 
field p neglected here. For additional sub-harmonic driv- 
ing the model predicts steady moving interfaces with the 
direction of motion controlled by the phase of the sub- 
harmonic component. Experimental study of the con- 
trolled interface motion will be reported elsewhere [pi . 
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FIG. 1. Phase diagram for Eq. (Q) at b = 4. 
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FIG. 2. (a-d) Interface instability and labyrinth formation, 
Eq. (0), uj = 2, b = 4,7 = 2.9, domain size 100 x 100 units, 
the snapshot are taken at times t = 1000, 1600, 2300, 4640; 
(e-h) labyrinth formation from a circular spot, SHE, 5 = 1.4, 
domain size 100 x 100 units, t = 200, 1300, 1700, 1900. 




FIG. 3. Interface velocity v for u) = 1, b = 4, 7 = 2.5 vs q 
at $ = 0. Inset: v vs $ at q = 0.01. (•) - numerical results, 
( — ) - analytical expression (|l7|) . 
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